Scaling behaviour in probabilistic neuronal cellular automata 
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We study a neural network model of interacting stochastic discrete two-state cellular automata 
on a regular lattice. The system is externally tuned to a critical point which varies with the degree 
of stochasticity (or the effective temperature). There are avalanches of neuronal activity, namely 
spatially and temporally contiguous sites of activity; a detailed numerical study of these activity 
avalanches is presented, and single, joint and marginal probability distributions are computed. At 
the critical point, we find that the scaling exponents for the variables are in good agreement with a 
mean-field theory. 

PACS numbers: 



I. INTRODUCTION 

Networks in a variety of natural settings are currently 
the subject of intense enquiry and span a variety of dis- 
ciplines, ranging from biology [IH3L. ecology [4[, engi- 
neering [H-0] to social interactions Q. The manner in 
which such networks arise, given the diversity of settings 
in which these are observed, is an additional issue. The 
idea that such structures are in some sense inevitable 
once a set of simple rules and interactions is specified, 
namely that there is "self-organization" at work has been 
a theme that has also been studied in diverse settings @ . 

The neural network structure — the nature and signif- 
icance of the interactions and connections in neuronal 
systems — has long been recognized and has been an area 
of study for several decades. There has been extensive 
work on the dynamical properties of neural networks, in 
particular those that emerge as a result of evolving con- 
nection strengths [ltj [H| , the response to external stim- 
uli [l2| , and the interplay between topology and intrinsic 
nodal dynamics [H|,[l4[ . A specific question is on how the 
response of a sensory system to external stimuli is opti- 
mized [HI, fl5T- fl7| . p articularly in the context of informa- 
tion processing [18 1 and whether the self-organized net- 
work is essential to this feature. When a neural network 
is in a critical state, the avalanches of spiking activity [lj| 
are critical, namely they are characterized by power-law 
distributions, and the role of network topology [2(| ap- 
pears to be limited. Studies of avalanche dynamics have 
been both theoretical [211-124} and experimental: for in- 
stance, it has been suggested that the propagation of 
spontaneous activity recorded from the rat cortex can be 
described in terms of critical avalanches with an event- 
size exponent of 3/2 (25l. l26j. 

How does the global dynamics of such networks change 
when the nature of the interaction is altered? This ques- 
tion is explored here through the use of both analytic 
and numerical tools and is the main focus of the present 
work. Similar issues have been examined in other recent 
works [27], HH although the geometry of the model and 
the dynamics differ significantly from the present case. 



In the following section, we describe the system studied 
here. This builds upon the model introduced earlier by 
Kinouchi and Copelli (KC) [l2[ by introducing stochas- 
ticity in the deexcitation process: excited neurons return 
to the silent state with some probability, reflecting the in- 
herently stochastic nature of cellular and subcellular dy- 
namics. In concordance with earlier work [12} there is a 
non-equilibrium phase transition, and here we determine 
the critical dynamics by computing the Lyapunov expo- 
nent. These results arc presented in Section III, and this 
is followed by a study of the avalanches of neuronal activ- 
ity in Section IV. A mean-field theory for the stochastic 
transition and a general methodology relating the scaling 
exponents are also given. The final Section V summarizes 
our results. 



II. THE GENERALIZED KC MODEL 

The neuronal cellular automaton model introduced by 
Kinouchi and Copelli is a quenched Erdos-Renyi network 
of excitable elements [Hj] . Nodes in the network can be in 
one of m states: (silent), 1 (excited or active) or to — 2 
refractory states (which are denoted — (m — 2), . . . , — 1 for 
convenience). The dynamics results from the following 
three processes: 

1. External driving or noise as a stimulus can cause a 
silent node to become active. This 0— >1 transition 
occurs at each time step with probability r\ = 1 — 
exp(— r), where r is the stimulus rate. 

2. At each time step, a silent node % is excited (namely 
— >1) if it receives a stimulus from one of its active 
neighbors j. This can occur with probability Aij, 
the weight of the connection. 

3. Spontaneous transitions that can be of two kinds, 
namely, an excited state becomes refractory ( 1— > 
— (to— 2) ) with unit probability and further, a node 
in state I moves to I + 1 if — (m — 2) < I < — 1, also 
with unit probability. 
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Network response F is calculated by the density of ac- 
tive sites pt averaged over a time window T. Measure of 
fluctuations in F is the susceptibility \ = (F 2 ) ~ (^) 2 j 
while the spontaneous activity at ?y = is denoted by 
Fq- This nonequilibrium system undergoes a continu- 
ous phase transition from a moving phase (correspond- 
ing to sustained neuronal activity) to an absorbing phase 
(where all nodes are in silent state) as a function of the 
average branching ratio (l2j . For every node in the net- 
work, the branching ratio gives the average number of 
excitations created by that clement at the subsequent 
time step. Similar interactions have been used to model 
epidemics on contact networks [2!| . 

The relationship between criticality and network struc- 
ture is crucial, since the response of this network to ex- 
ternal stimuli — and hence the dynamical range — is op- 
timal when the spontaneous activity of the network is 
critical |12[. However, in a generalization of the KC 



model that accounts for the integration of excitatory 
inputs [l7| . the order of the phase transition itself ap- 
pears to change: the transition becomes discontinuous 
and results in Instability in the system, with a history- 
dependent dynamical range apparently being optimized 
in this bistable regime. The largest eigenvalue A of the 
network connectivity matrix is the "universal" control 
parameter [20I ]. Furthermore, if the connectivity matrix 
is symmetric with nonzero elements all equal (namely all 
connections have equal weights) then A is proportional 
to the average degree of the network [3(| • 

We study a variation on the KC model on a two- 
dimensional square lattice using both the nearest neigh- 
bor (namely Von Neumann) as well as the Moore con- 
nectivity. For simplicity we take to = 2; nodes can be 
either silent or excited. Further, the external stimulus 
is absent, namely r\ = 0. Based on the electrophysio- 
logical differences, neurons can be broadly categorised as 
bursting or spiking. Considering a neuron as a discrete 
dynamical unit, excited state for a single time step char- 
acterizes the spiking behaviour whereas multiple consec- 
utive spikes of a bursting neuron preceding a refractory 
period are mimicked by a prolonged active state lasting 
several time steps [H, [l4j]. To account for both these 
behaviours in the present network, the transition from 
the active to the silent state (cf. process 3 in the above 
dynamics) is considered as a stochastic process. This 
transition is taken to occur with probability fi £ [0, 1], fi 
= being the "zero-temperature" limit, while fi — 1 cor- 
responding to the case studied earlier by KC for a spiking 
neuron. The elements of the symmetric connectivity ma- 
trix, the Ay's, are drawn from an uniform distribution 
in [0,1], and we rescale the largest eigenvalue A of this 
matrix appropriately to use as a control parameter poj . 
This rescaling is done in a way that the largest eigen- 
value of the new adjacency matrix lies between and 2 
with A = 1 being the critical point. The elements of this 
matrix are kept fixed throughout the simulation imply- 
ing quenched disorder. We choose to use A as the control 



parameter to further extend our work to the study of real 
networks with heterogeneous degree distributions where 
the branching ratio fails to predict criticality [3l|. Also, 
with the use of A as the control parameter, our calcula- 
tions become independent of the system size. 



III. RESULTS 

In the limit fj, — > 1, the present model reduces to that 
studied by KC. In order to further characterize the crit- 
ical dynamics we compute the Boolean Lyapunov expo- 
nent 



32j A by coevolving two initial network configura- 



tions and examining the average variation of the Ham- 
ming distance as a function of time. The Boolean Lya- 
punov exponent is given as 



1 n 
A=^Vlog 
77, z - 



(1) 



where n is the number of sample (random) perturbations 
considered, Tie is the number of time steps, and di n and 
d ou t are the initial and final Hamming distances. Suf- 
ficiently large Tie and n need to be taken for reliable 
results. Other definitions of the Boolean Lyapunov ex- 
ponents have also been extensively used earlier to charac- 
terise the dynamics of cellular-automata models [HI, [HJ . 

Our results are consistent with a critical point at A = 1 
for n = 1; nearby trajectories converged for sub-critical 
networks, namely A < for A < 1 while the dynamics is 
chaotic in the super-critical networks: A > for A > 1. 
Shown in Fig.QTa) are the numerical results at this value 
of /i showing that both A and A-l cross the axis at the 
same point. 

For however, we find that the critical point is 

no longer at A = 1 and is instead shifted by an amount 
proportional to fi. The change in the sign of the Boolean 
Lyapunov exponent, the phase transition in Fq, and the 
optimization of the susceptibility shown in Fig. [lja)-(c) 
respectively; all three indicators clearly show a shift in 
the critical point for different values of /i. Similar re- 
sults are obtained for m > 2; wherein a node in state 
1 (excited) goes to state —(to — 2) (refractory) with a 
probability [i and its subsequent updates till it reaches 
the state (silent) are deterministic; an extended per iod 
of refractoriness does not affect the critical point [35| in 
the present setting. In ID (namely K = 2), the system 
does not show a phase transition at /i = 1 for a 3 state 
model (to = 3) but for [i < 1, the phase transition is 
recovered [36j . 

We confirm the linear variation of the critical value, 
namely A c ~ y, both numerically and analytically us- 
ing mean-field analysis. In the KC model the transition 
from the silent to the active state (namely — > 1) can 
take place in two ways: (i) through an external stim- 
ulus, and (ii) via active neighbor interactions. On a 
two dimensional square lattice, the degree distribution is 
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FIG. 1: (a) Boolean Lyapunov exponent A for Tie = 10 4 and n — 10 3 , (b) Spontaneous activity Fo averaged for 10 3 time steps 
and (c) Susceptibility \ with control parameter A for different values of /i for 10 6 nodes, 77 = and Moore neighborhood. 





FIG. 2: Plot of avalanche size distribution P(S) for different 
A for 10 7 avalanche events. 



FIG. 3: Joint size-duration probability distribution P(T, S) 
for total nodes, iV = 1024 2 using log binned data. 



P(k) = 8(k — K), where K denotes degree of each node. 
The probability of causing a silent node to become active 
at any time t via the input it receives from at least one of 
its active neighbors is given by pt = 1 — [1 — (\/K)F t ] K 
where F t denotes the fraction of active nodes at time t 
and in the mean field approximation (Aij) is X/K. 

The transition from state 1 — > occurs with probability 
fi. If Pi (t) denotes the configuration probability of active 
nodes and Po(t) that of silent nodes, the corresponding 
master equation is 

Pi(t + 1) = T}Pb(t) + (1 - v)PtPo{t) + (1 - (2) 

Denoting Pi (t) = F t , one has Po (i) = 1 — F t and in terms 
of Ft, the master equation becomes 

F t+ i = [1 - F t ][r, + (1 - V )Pt] + (1 - n)F t . (3) 
In the stationary state F t +i = F t = F, giving 

A 



»F = (1 - F) 



1 - (1 - 77) 1 



K 



:F 



K 



(4) 



In the absence of an external stimulus namely 77 = 0, 
by setting A as A c + e, in the limit e —> 0, F —> we get 



A c -> n, 



and 



F{\) 



A — /i 

C ' 



(5) 



(6) 



where the constant C = fx + ^IT^ 2 - Therefore, since 
F{\) scales as |A — \ c \@ , we get /3 = 1, although the 
numerical results (f3 num = 0.499 ± 0.002) are not in good 
agreement with this estimate. 

At the critical point A c = \i and with external stimulus 
77, the resulting behaviour is 



F(r) a v^/C 



(7) 



where C is an arbitrary constant. The numerical results 
are in better agreement with this behaviour, namely the 
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FIG. 4: The data collapse of Pt(S) curves for the joint size- 
duration distribution. Each Pt (S) curve is the intersection of 
the joint size-duration distribution P(T, S) with the T — To 
plane. The collapse width is w = 4.81 xlCP 2 . 



scaling F(r) ~ r 1 / 5 with 6 = 2 « 5 n «m (not shown here). 
In the limit A — /i — > e, near the critical point, _F — > 0, 
wc find that the susceptibility varies as 



(8) 



namely the scaling \ ~ l e | 7 ( £ > 0) an d X ~ l e l 7 ( e < 
0) with 7 = 7' = 1. Numerically, however, we find 7„ um 
= 0.89 ± 0.04, i num = 1.15 ± 0.04. 

A pertinent observable that has been introduced for 
such excitable systems is the dynamic range 



A = 10 log! 



rp.g 
ro.i 



(9) 



where the external stimulus r x leads to response F x . A 
measures the range of stimuli that can be discerned by 
the changes in the network response F. Previous studies 
have reported that the dynamic range attains its max- 
imum value A max at the critical point of the system 
[H |!|. With the addition of stochasticity, A max oc- 
curs at the shifted critical point consistent with Fig. [TJ 
however, its magnitude appears to be independent of /i 
(cf. Eqs. © and ©). 

The mean-field exponents obtained so far are /3 = 1, 7 
= 1, and 6 = 2. Now, we consider the cluster size prob- 
ability distribution which obeys P(S) ~ S^ r+1 - S~ TS 
(for details refer section ITV|) . Below the critical point , 
the cutoff cluster size is related to correlation length 



cLS Sq 



-l/a 



£ , where e is the con- 



trol parameter and D is the fractal dimension of cluster. 
The order parameter F is characterized by 



dSP(S) 



(10) 



where P(S) ~ 5 _r+1 /(5'|e| 1 / cr ) and / a general scaling 
function. Changing variables to y = S'|e| 1 /' T , the above 
integral reduces to 



|(r-2)/<x 



dyy- T+L f(y), 



(11) 



and hence on comparison with F ~ \s\P, we infer that 
/? = (t — 2)/er. The different moments of P(S) connect 
the various exponents via simple algebraic equations, and 
with a little computation, it can be seen that a = 3 — 2/ a, 
P =(r-2)/a, 7 = (3-r)/a, D = 1/W, and r = 2 + 1/6 
[13, HI]- Using the known values of 6 and 7, it is a simple 
matter to get r = 5/2 , a = 1/2, and a = -1. 

If we assume that the avalanche dynamics on the net- 
work is diffusion like [39| , then since the size of the largest 
cluster is limited by the linear system size L, any pertur- 
bation at the critical point takes an average L 2 steps to 
diffuse away. So, from the stationarity condition we get, 
7/1/ = 2, resulting in v = 1/2 [3^| and D = 4. Using the 
hyperscaling relation D = d — [3/v wc get d = 6 [4(| as 
the upper critical dimension. 



IV. 



AVALANCHE PROPERTIES 



Starting with a null configuration (all nodes are in 
state 0) we select a site at random and make it active. 
This random flipping induces activity in the neighbor- 
hood which may then propagate, but this will eventually 
die out so that the system returns to the null configura- 
tion. This entire event constitutes an avalanche and the 
total activity at each time step is the avalanche signal 
s(t). If T is the total time, then the total activity during 
this event 



S = $>(i), 



(12) 



is the avalanche size S. We also compute other re- 
lated quantities, the area A which is equal to the total 
number of distinct sites which become active during the 
avalanche, the energy 



E 



and the magnitude 

M = max|s(f)| : < t < T. 



(13) 



(14) 



The probability distributions of each of these quanti- 
ties have been computed; shown in Fig.[5]is the avalanche 
size distribution P(S) for different values of A where the 
power-law behaviour in the super-critical, critical and 
sub-critical regimes can be clearly identified. The scaling 
relations for these distributions can be defined in general 
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as P(X) <~ X~ Tx , where X can be size 5, duration T, 
area A, magnitude M or energy E. 

The various exponents are, of course, related. If X, Y 
are a pair of the above variables and if (X) ~ (Y)~ (XY , 
where (•} denotes the average value, then there is a simple 
relation among the exponents tx, ty and ~/xy, 



7XY 

Further, using the relationship (Y) ~ X lYX , one gets 

T X = 1 - (1 - Ty)7yx, (16) 

and from Eq. (JTSJ) and (fPoj) . 

7xr = — ^— • (17) 

Tyx 

The scaling form for the joint probability distribution, 
P(X, Y), can be obtained via 

P(X ax X,X aY Y) = XP(X,Y). (18) 

where, ax and ay are the corresponding exponents asso- 
ciated with the arguments of a generalized homogeneous 
function. 

By putting A ~ X~ 1 / ax in Eq. ([18]). we get 

P(X,Y) ~ X^P(Y/X aY/ax ). (19) 
Comparing this to the expression 

P(X,Y) ~ X-( TX ^ X )P(Y/X' ryjr ), (20) 
one can see that 

ay 1 

7rx = , = -(tx + 7yx)- (21) 

a x ax 

The replacment X — > Y and Y -> X in Eq. (|2T|) simi- 
larly gives 

1 

7xy = — , — = -(ry + 1XY)- (22) 
ay ay 

Clearly, Eqs. (f2"T|) and (|22|) can be used to obtain Eqs. ([Toi 

We study networks with degree K = 4 and system 
size iV = 1024 2 at the critical point (for fi = 1) and 
compute joint distributions |4ll l42l | for four pairs of the 
above related quantities, size-duration (see Fig. [3]) as well 
as size-magnitude, size-area and size-energy (which are 
not shown here). For each of these joint probability dis- 
tributions P(X, Y) we fix one variable Y to obtain the 
marginal- X distribution, Py(X). This obeys the scaling 
relation Eq. (|20|) . 

A typical data collapse of the marginal- A~ distribu- 
tions is shown in Fig. 01 and critical exponents obtained 
through this analysis are given in Table HI As can be seen, 



Exponent Measured value Scaling relation 



TS 


1.49 ± 0.03 




TT 


1.82 ±0.06 


1.77 


TA 


1.69 ±0.06 


1.72 


TM 


1.97 ±0.07 


2.10 


TE 


1.35 ±0.02 


1.35 


7ST 


1.527 ±0.002 




■jsa 


1.474 ± 0.003 




7SM 


2.25 ±0.03 




7SB 


0.7712 ± 0.0005 





TABLE I: Exponents characterizing avalanche properties. 
Measured values are obtained from data collapse while the 
values mentioned in the third column are obtained from the 
scaling relation Eq. (|15[) . 



these are in good agreement with the corresponding ex- 
ponents obtained through the scaling relations (within 
statistical error). The collapsing width for the vari- 
ous curves was also calculated [41(. Although in the 
present study the network is externally driven to crit- 
icality by varying the control parameter A, avalanches 
showing power-law behaviour also form an integral part 
of the study of self-organized critical neuronal systems 
0- 



V. DISCUSSION AND SUMMARY 

Recent efforts to understand the operation of neuronal 
systems have examined a number of models wherein the 
feature of self-organization is important. This would be 
one way of understanding the emergence of complexity in 
the context of neuroscience. Although automaton mod- 
els are by their very nature caricatures of how a biolog- 
ical neuron operates, they capture some of the essential 
features. Furthermore, the manner in which neuronal 
communication operates can also be modeled within the 
framework of coupled automata. 

Here we investigated the automaton model introduced 
by Kinouchi and Copelli [l2j ■ While it is known that the 
critical point of this model is invariant under the change 
of topology [n| , we alter the nodal dynamics of the KC 
model and consider its impact on criticality. Other recent 
studies have also discussed the effect of these variations 
in a different context (28j . underscoring the considerable 
current interest in such problems. We incorporate ad- 
ditional stochasticity in the deexcitation transition from 
the active to the silent state and the primary effect of 
this is to shift the point of criticality without destroying 
it. This extended model thus retains the main features 
of criticality and emergence. Although the effect of hav- 
ing a variable number of refractory states on the critical 
point in the KC model has been investigated [35| , in the 
present case, the dynamics tends to hold the neuron in 
the excited state for a longer duration thus taking into 
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account the bursting nature as well. Furthermore, ow- 
ing to the inherent stochasticity in the intrinsic dynam- 
ics of a neuron, the duration of bursts may also differ 
([43| and references therein). Indeed, experiments sug- 
ges t that there exists variability in the burst durations 

Activity in the network is characterized through 
avalanches, namely sets of firing nodes that arc 
connected spatially and temporally. Probability dis- 
tributions for a variety of quantities that quantify the 
avalanches can be described through a set of exponents 
which we have computed in a two dimensional realization 
of this model. We also developed a mean-field model 
for the dynamics, and show that at the critical point 
the numerical results are in good agreement with the 
theoretical predictions, and further, with the exponents 
obtained for branching processes [45|, |4(| . The KC model 



established a theoretical framework for the optimization 
of dynamical range at criticality, and this has been 
subsequently verified experimentally [47l l48j . It would 
also be of interest to see if optimization of susceptibility 
can also be demonstrated experimentally in such models. 
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